﻿using LxBasic;
using LxLESS;
using System;
using System.Diagnostics;

namespace LxCore
{
    public class MassFluxBase : ZoneHead, ICalResidual
    {
        //local
        public double[] MassImbalance;//向外流出为正
        public double[] InnMassFlux;
        public double[] BouMassFlux;
        protected bool isMassFluxInitialized { get; set; } = false;
        protected bool isMassImbalanceInitialized { get; set; } = false;
        //share
        public IField<double> Source = ConstField.Double0;
        public override void CheckCreate(FvZone zone)
        {
            base.CheckCreate(zone);
            bool newInnMf = CheckCreator.ArrayStrict(ref InnMassFlux, InnFaces.Length);
            bool newBouMf = CheckCreator.ArrayStrict(ref BouMassFlux, BouFaces.Length);
            if (newInnMf || newBouMf)
            {
                isMassFluxInitialized = false;
                isMassImbalanceInitialized = false;
            }
            bool newImb = CheckCreator.ArrayStrict(ref MassImbalance, InnCells.Length);
            if (newImb)
            {
                isMassImbalanceInitialized = false;
            }
        }
        public override void Reuse(ZoneHead oldOne)
        {
            base.Reuse(oldOne);
            if (oldOne is MassFluxBase o)
            {
                this.MassImbalance = o.MassImbalance;
                this.InnMassFlux = o.InnMassFlux;
                this.BouMassFlux = o.BouMassFlux;
            }
        }
        public override void Initialize()
        {
            if (!isMassFluxInitialized)//TODO 如果是读取文件获得了Velocity，需要计算MassFlux
            {
                //MultiThreads.RunParallelTask((threadInfo) =>
                //{
                //    LxTool.ClearArray(InnMassFlux, threadInfo);
                //    LxTool.ClearArray(BouMassFlux, threadInfo);
                //    LxTool.ClearArray(MassImbalance, threadInfo);
                //});
                //isMassFluxInitialized = true;
                //isMassImbalanceInitialized = true;
                MultiThreads.RunParallelTask(CalFaceMassFlux);
                isMassFluxInitialized = true;
            }
            if (!isMassImbalanceInitialized)//如果读取的文件，此时没有massImb数据
            {
                MultiThreads.RunParallelTask(CalMassImbalance);
                isMassImbalanceInitialized = true;
            }
        }
        public virtual void CalFaceMassFlux(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
        }
        public virtual void CalMassImbalance(ThreadInfo threadInfo)
        {
            (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            if (Source.IsEqual0())
            {
                for (; i1 < endp1; i1++)
                {
                    MassImbalance[i1] = InnCells[i1].CalIfOwnerAddElseMunus(InnMassFlux, BouMassFlux);
                }
            }
            else
            {
                var s = Source.GetCellVarArr(Zone.Inner);
                for (; i1 < endp1; i1++)
                {
                    MassImbalance[i1] = InnCells[i1].CalIfOwnerAddElseMunus(InnMassFlux, BouMassFlux) - s.GetCellVar(i1);
                }
            }
        }

        #region ICalResidual
        private Float128 SumOfImbRateSquare = 0;
        public unsafe void CalResidual(ThreadInfo threadInfo)
        {
            {
                (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
                Float128 sum_Thread = 0;
                for (; i1 < endp1; i1++)
                {
                    double rate = MassImbalance[i1] / (calMassFlux(ref InnCells[i1]) + 1e-100);
                    sum_Thread += rate * rate;
                }
                lock (MultiThreads.GlobalLock)
                {
                    SumOfImbRateSquare += sum_Thread;
                }
            }
            double calMassFlux(ref FvInnCell c)
            {
                double flux = 0;
                fixed (int* start = c.AdjFaces)
                {
                    int* adjf = start, endp1 = start + c.InnNFaceEnd;
                    for (; adjf < endp1; ++adjf)
                    {
                        flux += Math.Abs(InnMassFlux[*adjf]);
                    }
                    endp1 = start + c.AdjFaces.Length;
                    for (; adjf < endp1; ++adjf)
                    {
                        flux += Math.Abs(BouMassFlux[*adjf]);
                    }
                }
                return flux;
            }
        }
        public (bool isValid, double res) GetResidualReset()
        {
            double res = Math.Sqrt((double)SumOfImbRateSquare / InnCells.Length);
            SumOfImbRateSquare = 0;
            return (true, res);
        }

        #endregion
        public override void FromJsonKV(string key, IJsonReader jsReader)
        {
            switch (key)
            {
                case nameof(InnMassFlux): jsReader.ReadDoubleArray(InnMassFlux, e => InnMassFlux = e); break;
                case nameof(BouMassFlux):
                    jsReader.ReadDoubleArray(BouMassFlux, e => BouMassFlux = e);
                    isMassFluxInitialized = true;
                    break;
                default: base.FromJsonKV(key, jsReader); break;
            }
        }
        public override void ToJsonKV(IJsonWriter jsWriter)
        {
            base.ToJsonKV(jsWriter);
            //需要存储massFlux，massFlux的计算用到SIMPLE Pressure中的Dc，而Dc又受到Velocity方程系数的影响，方程系数又受到massFlux影响
            if (isMassFluxInitialized)
            {
                jsWriter.WriteDoubleArray(nameof(InnMassFlux), InnMassFlux);
                jsWriter.WriteDoubleArray(nameof(BouMassFlux), BouMassFlux);
            }
        }
    }
}
//Array.Clear(InnVelDirCor, i1, endp1 - i1);
//Array.Fill(InnVelDirCor, 1, i1, endp1 - i1);
//protected float calVelDirCor(double density, Vector3d v_O, Vector3d v_N, double OToNLength)
//{
//    double u = Vector3d.DotProduct(v_O, v_N) > 0 ? Vector3d.CrossProduct(v_O, v_N).Length() : Math.Sqrt(v_O.LengthSquared() * v_N.LengthSquared());
//    return (float)(density * u * OToNLength);
//    //return (float)(density * OToNLength * (v_O - v_N).Length() *0.5* (1 - Vector3d.DotProduct(v_O, v_N) / (Math.Sqrt(v_O.LengthSquared() + v_N.LengthSquared()) + 1e-100)));
//}
//protected float calVelDirCor(Vector3d v_O, Vector3d v_N)
//{
//    return (float)Math.Max(0, Vector3d.DotProduct(v_O, v_N) / (Math.Sqrt(v_O.LengthSquared() * v_N.LengthSquared()) + 1e-100) * 10 - 9);
//}